#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

sig.effects(measure = "FA", atlas = "glasser", dataset = "pnc")

sig.effects(measure = "FA", atlas = "glasser", dataset = "hcpd")

#Function to calculate the number and percent of thalamocortical connections showing a significant age*sex interaction
sig.effects.sex <- function(measure, atlas, dataset){
  sexeffects.df <- get(sprintf("sexeffects_%s_%s_%s", measure, atlas, dataset)) 
  sexeffects.df$significant <- p.adjust(sexeffects.df$GAM.int.pvalue, method = c("fdr")) < 0.05 #fdr-corrected 
  sigeffect.totaln <- sum(sexeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(sexeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age by sex interactions", dataset, sigeffect.totaln, sigeffect.percent))
}

sig.effects.sex(measure = "FA", atlas = "glasser", dataset = "pnc")

sig.effects.sex(measure = "FA", atlas = "glasser", dataset = "hcpd")

smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)

smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")

ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  scale_y_continuous(limits = c(-0.02, 0.0138)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)

#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
  
  #Read in participant-level data and model-level fitted values and derivatives 
  participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
  fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA

  #Connection spline plot with participant-level data
  trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
    geom_point(color = alpha(display_color, 0.5), size = .1) +
    geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
    geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted,  ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
    labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
    theme_classic() +
    scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    scale_y_continuous(breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  #Significant derivatives (rate of age-related change) plot 
  derivatives.plot <- ggplot(data = derivatives.df) + 
    geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
    scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color,  na.value = "white") +
    labs(x="\nAge", fill = "Derivative") + 
    theme_classic() +
    scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    theme(
        legend.position = "none",
        axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size = 6, family = "Arial", color = c("black"))) 
  
  allplots <- list(trajectory.plot, derivatives.plot)
  tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)

  return(tract.plot) 
}

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") 

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.06, 0.15)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") 

ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
  geom_point(color = c("#9c3a80"), size = 0.4) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
  scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Agematuration_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.94, height = 1.58)

cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)

cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late

matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey

ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Maturationalbrackets_PNC.pdf", device = "pdf", dpi = 500, width = 4.4, height = 4.4)

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))

devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.pnc

gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")

devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
  arrange(desc(term.correlation))

devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))

devpattern.neurosynth.hcpd

merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)

ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 2.2, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
  scale_y_discrete(
  labels = c("expectancy"=expression(bold(expectancy)),
  "monitoring"=expression(bold(monitoring)),
  "cognitive_control"=expression(bold(cognitive_control)),
  "reasoning"=expression(bold(reasoning)),
  "response_inhibition"=expression(bold(response_inhibition)),
  "thought"=expression(bold(thought)),
  "recall"=expression(bold(recall)),
  "induction"=expression(bold(induction)),
  "object_recognition"=expression(bold(object_recognition)),
  "visual_perception"=expression(bold(visual_perception)),
  "gaze"=expression(bold(gaze)),
  "induction"=expression(bold(induction), parse=TRUE))) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_line(color = c("gray90")),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)

#Create 10,000 spins of the HCPD age of maturation map for a neurosynth term overlap null analysis. Each spatially rotated map will be correlated with neurosynth term maps as above and the 10 most positively and negatively correlated terms will be identified. The number of overlapping terms with the empirical PNC terms will then be computed, creating a term overlap null based on spun maps for computing a pspin value.
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd) #complete data for spinning
x <- spin.data$smooth.increase.offset #specific cortical map to spin

#spin the empirical data 
perm.id <- perm.id.full
nroi = dim(perm.id)[1]  #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) { #spin it
  for (i in 1:nroi) {
    x.perm[i,r] = x[perm.id[i,r]]
  }
}
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm$label <- spin.data$label

#compute neurosynth correlations for all spins and identify overlapping terms
compute.termoverlap <- function(perm){
  this.perm <- sprintf("perm%s", perm)
  nullmap.data <- x.spatialperm %>% select(label, all_of(this.perm))
  nullmap.neurosynth <- left_join(nullmap.data, neurosynth.terms, by = "label")
  nullmap.neurosynth <- left_join(nullmap.neurosynth, glasser.regions, by = "label")
  
  #null neurosynth correlations for this spin
  devpattern.neurosynth.null <- ldply(neurosynth.termlist, function(t){
    developmentmap.neurosynth.decoding(df = nullmap.neurosynth, dev.metric = this.perm, term = t)}) %>% 
    arrange(desc(term.correlation))
  
  #identify top and bottom most correlated terms for this spin
  devpattern.neurosynth.null <- rbind(slice_max(devpattern.neurosynth.null, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.null, order_by = term.correlation, n = 10))
  
  #compute overlap
  pnc.terms.empirical <- devpattern.neurosynth.pnc %>% select(term)
  null.terms <- devpattern.neurosynth.null %>% select(term)
  num.overlap <- length(intersect(pnc.terms$term, null.terms$term))
  
  return(num.overlap)
}

term.overlap.null <- lapply(1:nperm, compute.termoverlap)
term.overlap.null <- do.call(rbind, term.overlap.null) %>% as.data.frame %>% set_names("overlapping.terms")

write_rds(term.overlap.null, "/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/termoverlap_null.rds")

term.overlap.null <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/termoverlap_null.rds")

sprintf("In %s out of 10,000 spins, the term overlap number was equal to or greater than 11", sum(term.overlap.null$overlapping.terms >= 11))

#calculate spin based p
sum(term.overlap.null$overlapping.terms >= 11) / 10000

derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)

derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") #fdr-correct this for axis-based analyses below

ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey

ggplot() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  17.5, na.value = "gray93", limits = c(16, 19), oob = squish) + 
  new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) 

S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)

ggplot() + 
  geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  180, na.value = "gray93") + 
   new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) 

cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") #fdr-correct this for axis-based analyses below

ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_hcpd.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

#Glasser parcel x,y,z coordinates for defining A-P, M-L, and D-V axes
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$x.mni[1:180]  <- hcp_mmp1.0$x.mni[1:180]*-1 #convert R-> L hemisphere coords into medial lateral coords
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
hcp_mmp1.0 <- hcp_mmp1.0 %>% select(-hemi)

#C-Mt files for defining the group-level core-matrix gradient
CPt.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv") %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv")  %>% select("tract", "label", "orig_parcelname", "mean_CPt")
CPt.gradient <- merge(CPt.glasser.pnc, CPt.glasser.hcpd, by = c("label", "tract", "orig_parcelname"), all.x = TRUE, all.y = TRUE, suffixes = c("_pnc", "_hcpd"))
CPt.gradient$CPt <- CPt.gradient %>% select(mean_CPt_pnc, mean_CPt_hcpd) %>% rowMeans(na.rm = T)
CPt.gradient <- CPt.gradient %>% select(CPt, label)

hcp_mmp1.0 <- left_join(hcp_mmp1.0, CPt.gradient, by = "label")

#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
  df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
  df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes,  by = c("label", "tract", "orig_parcelname"))
  
  #S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - development correlation
  anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "CPt")

axis.cor.ps <- c(0.00075, 0.0621, 0.43765, 0.26005, 0.0287) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))

cocor.ps <- c(0.0001552829, 3.832011e-08, 0.0002034701, 0.0006272672)
p.adjust(cocor.ps, method = c("fdr"))

#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4877183)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3101209)
axis.results <- axis.results %>% add_row(Axis = "C-M", corr = .2897652)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1173292)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = 0.03024772)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "C-M", "D-V", "M-L"))

ggplot(axis.results, aes(x = Axis, y = corr)) + 
  geom_segment(aes(x = Axis, xend = Axis,  yend =  
                     corr, y = 0), linewidth = 0.2) +
  geom_point(aes(fill = Axis), size = 3.6, shape = 22, color = alpha("white", 0)) +
  scale_y_continuous(limits = c(-0.01, .5)) +
  labs(x="Neuroaxis") +
  labs(y="\nAxis Alignment") +
  theme_classic() +
  scale_fill_manual(values=c("#672975", "#9898BF", "#9898BF", "#9898BF", "#9898BF")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.text.x = element_text(vjust = 0.7, angle = 20),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_lollyplot.pdf", device = "pdf", dpi = 500, width = 1.36, height = 1.05)

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "CPt")

axis.cor.ps <- c(0.00155, 0.02195, 0.3695, 0.38565, 0.0163) #S-A, A-P, M-L, D-V, C-M
p.adjust(axis.cor.ps, method = c("fdr"))

cocor.ps <- c(0.124339, 5.601018e-07, 1.390815e-05, 0.02169413)
p.adjust(cocor.ps, method = c("fdr"))

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)

#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)

#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full

library(ciftiTools)

#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels

#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version 
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released 
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)

#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")

#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)

detach("package:ciftiTools", unload=TRUE)

EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)

df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

plasticity.dev$agedecline.protracted <- plasticity.dev$amplitude.agedecline
plasticity.dev$agedecline.protracted[plasticity.dev$agedecline.protracted == "NaN"] <- 23 #lets plot regions with most protracted functional expression of plasticity (amplitude never declines!)
snr.exclude <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/SNRmask_glasser360.csv") %>% filter(SNR.mask == 0) #but not those we excluded for low SNR
plasticity.dev$agedecline.protracted[plasticity.dev$label %in% snr.exclude$label] <- NA

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = agedecline.protracted, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(10, 17), oob = squish, na.value = "white") 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.7, height = 5.4) 

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 1), alpha("#323280", 0.8), alpha("#323280", 0.6), alpha("#323280", 0.4), alpha("#323280", 0.2)), limits = c(0.012, 0.016), oob = squish, na.value = "white") 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.4, height = 5.13) 

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c(alpha("#323280", 0.2), alpha("#323280", 0.4), alpha("#323280", 0.6), alpha("#323280", 0.9), alpha("#323280", 1)), limits = c(-0.04, -0.035), oob = squish, na.value = "white") 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 5, height = 5.08) 

plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)

EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)

myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)

amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

plasticity.ps <- c(EI.thalamus.pnc, myelin.thalamus.pnc, amplitude.thalamus.pnc)
plasticity.ps

p.adjust(plasticity.ps, method = c("fdr"))

plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)

EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)

myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)

amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#A2A2C8") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(10, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

plasticity.ps <- c(EI.thalamus.hcpd, myelin.thalamus.hcpd, amplitude.thalamus.hcpd)
plasticity.ps

p.adjust(plasticity.ps, method = c("fdr"))